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ABSTRACT 

A direct solution procedure for computing the flutter Mach number and the flut- 
ter frequency is applied to the aeroelastic analysis of propfans. The proce- 
dure uses a finite element structural model and an unsteady aerodynamic model 
based on a three-dimensional, subsonic, compressible lifting-surface theory. 

An approximation to the Jacobian matrix that improves the efficiency of the 
iterative process is presented. The Jacobian matrix is indirectly approximated 
from approximate derivatives of the flutter matrix. Examples are used to 
illustrate the convergence properties. The direct solution procedure facili- 
tates the automated flutter analysis and contributes to the efficient use of 
computer time as well as the analyst's time. Further details of the numerical 
procedure are given by Murthy and Kaza (1987). 
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MOTIVATION 


Flutter of propfans and other types of turbomachinery blading is an important 
phenomenon that has generated considerable interest. Flutter prevention has 
been a significant factor in the design of propfan blades. Flutter prevention 
is also significant for turbomachinery, particularly for unshrouded blades. 
With the recent advances in computer technology, automated design of propfan 
and turbomachinery blades by using optimization techniques has become practi- 
cal. Design optimization employing flutter constraints requires repeated solu 
tion of the aeroelastic equations of motion to obtain the flutter parameters 
as the design is updated. For the optimization to be performed in a realistic 
period of time an automated flutter analysis capability is essential. It is 
also desirable for the analysis to be computationally efficient in order to 
keep the central processing unit (CPU) time and the turnaround time within rea 
sonable limits. Automated flutter analysis can also shorten the nonautomated 
design process by reducing the analyst's time. 


WHY AUTOMATED FLUTTER ANALYSIS? 

• ESSENTIAL FOR REPEATED EXECUTION 
OF FLUTTER ANALYSIS CODE 



• ALSO USEFUL IN NONAUTOMATED DESIGN PROCESS 
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FORMULATION OF FLUTTER ANALYSIS PROBLEM 

The computational procedure to be presented is applied to the analytical formu- 
lation described in detail by Kaza, et al. (1987a, 1987b). This formulation 
is applicable to the flutter analysis of a single-rotor propfan containing an 
arbitrary number of blades rotating at a fixed speed in an axial flow. The 
structure is modeled by finite elements. The aerodynamic model (Williams and 
Hwang, 1986) is based on a three-dimensional, subsonic, compressible 
lifting-surface theory. 

For simplicity, the effect of steady deformations due to aerodynamic loads on 
the flutter boundary is neglected. The error introduced by ignoring the steady 
aerodynamic deformations is shown in Kaza, et al. (1987a) to be small enough to 
produce an approximate flutter point. In particular, the effect of the steady 
aerodynamic deformations on the flutter Mach number is not very significant. 
Thus the approximate flutter analysis neglecting steady aerodynamic deforma- 
tions is suitable for use in design optimization procedures that require 
repeated efficient execution of the flutter analysis. The optimal design can 
be easily checked for the flutter condition by using the refined flutter analy- 
sis with steady deformations and the conventional procedure. 

The propfan is assumed to have identical groups of blades symmetrically distri- 
buted about a rigid disk. The linearized aeroelastic equations of motion are 
then uncoupled for different intergroup phase angle modes o r . The flutter 
Mach number for the propfan is then the lowest Mach number at which one of the 
intergroup phase angle modes becomes unstable. 


I 

i 


[Mg][q] + [K g ](q] = [A(M,co)][q] 

[M g ] GENERALIZED MASS MATRIX 

[K g ] GENERALIZED STIFFNESS MATRIX 

(q) GENERALIZED COORDINATE VECTOR 

[A(M,«)J GENERALIZED AERODYNAMIC MATRIX 


[A(M,co)] IS USUALLY VALID ONLY FORSIMPLE HARMONIC MOTION 
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CONVENTIONAL PROPFAN FLUTTER ANALYSIS 


The conventional procedure for obtaining the flutter Mach number Mp and fre- 
quency cop is as follows. The aerodynamic matrix is evaluated at an assumed 
Mach number Mp, an assumed frequency cop, and an assumed interblade phase 
angle ap, and then the flutter equation is solved for all the eigenvalues 
This procedure has been implemented for propfans in a program called ASTR0P3. 

In general, these eigenvalues are complex. The real and imaginary parts of 
ico represent the effective damping and frequency, respectively. The assumed 
frequency is varied until it is equal to the frequency corresponding to the 
eigenvalue with the least effective damping. This frequency matching forms the 
inner iteration. When this iteration reaches convergence, the Mach number is 
varied until the effective damping of the eigenvalue corresponding to the 
matched frequency is equal to zero. This forms the outer iteration. The flut- 
ter Mach number and the flutter frequency are obtained at the convergence of 
the outer iteration. 

The conventional procedure cannot be reliably automated because it requires 
that the identity or ordering of the eigenvalues be preserved over a wide range 
of assumed frequencies and Mach numbers. Most eigensolution routines do not 
compute the eigenvalues in any particular order, and the sorting of eigenvalues 
by frequency or magnitude does not usually preserve the continuity. Loss of 
continuity necessitates user intervention and complicates the automated analy- 
sis. A direct solution of the flutter equation that alleviates these problems 
is proposed and described. It views the flutter equation as an inplicit 
double-eigenvalue problem. 



• NEEDS EIGENVALUE TRACKING. AND IS THUS DIFFICULT TO AUTOMATE 

• REQUIRES DOUBLE ITERATION (INNER-OUTER LOOPS) 

• NEEDS AS MANY EIGENVALUES AS THERE ARE MODES 
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DIRECT SOLUTION OF FLUTTER EQUATION 


When the dependence of the aerodynamic matrix on the assumed Mach number and 
frequency is considered explicitly, the flutter equation can be viewed as an 
implicit double-eigenvalue problem. In general, the aerodynamic matrix [A(Mp, 
cop] is a transcendental function of the assumed frequency and Mach number. 

Only real values of Mp and cop are of interest. The two eigenvalues Mp 
and cop are coupled since the aerodynamic matrix is complex. 

We now have two equations in two unknowns, Mp and cop. These equations can 
be solved by any of the methods for solving a system of nonlinear equations. 
When they are solved for Mp and cop, no further iterations are required for 
the purpose of matching assumed and computed quantities. This procedure is 
illustrated below. Once Mp and cop are found, inverse iteration can be 
used to find the flutter mode. 


IF Mp IS THE FLUTTER MACH NUMBER AND cop THE FLUTTER FREQUENCY, 

Who) = (0) 


WHERE 


[B] = - Wp[Mg] + [Kg] - [A(Mp,U F )] 


FOR A NONTRIVIAL FLUTTER MODE, WE HAVE 


det[-<of[M g ] + [K g ] - [A(M f ,« f )]] = 0 

D = det(-cop[M g ] + [K g ] - [A(M F ,w F )]] 
= Op(M F ,co F ) + iD|(M F ,w F ) 


WHERE Dr AND D ( ARE THE REAL AND IMAGINARY PARTS OF THE CHARACTERISTIC 
DETERMINANT D, RESPECTIVELY. THEN AT FLUTTER CONDITION 


Dr(M F ,« F ) = 0 
D|(Mp,up) = 0 
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PROPFAN FLUTTER ANALYSIS BY DIRECT SOLUTION 


In contrast to the conventional procedure the direct solution outlined here 
eliminates the need to track eigenvalues to determine the flutter point. In 
addition to this important benefit the double iteration on the complex eigen- 
solution is replaced by a single solution of a system of two real nonlinear 
equations. The Mach number and frequency are varied simultaneously in this 
procedure rather than one at a time as in the conventional procedure. Thus the 
flutter Mach number and the flutter frequency are determined simultaneously. 

The formulation of a transcendental double-eigenvalue problem in preference to 
a linear single-eigenvalue problem may seem to defeat the objective of 
increased efficiency, even if it is more suitable for automation. However, the 
price to be paid is not as great as it may seem. The transcendental eigenvalue 
problem needs to be solved for only one set of eigenvalues in most cases, 
whereas the linear eigenvalue problem has to be repeatedly solved for all the 
eigenvalues, which are equal in number to the number of assumed mode shapes. 

The direct solution may not find the lowest flutter Mach number if more than 
one structural mode were to flutter in the Mach number and frequency range of 
interest for the selected intergroup phase angle mode. Under these circum- 
stances one will be forced to search the entire range of interest for the roots 
Mp and top, starting with different initial guesses. This is not amenable to 
an efficient automated procedure. However, it is expected that these circum- 
stances will rarely occur for tuned or alternately mistuned propfans. This is 
not a major limitation for two other reasons: (1) the frequency interval in 

which flutter occurs is usually determined early in the design phase and (2) 
the search domain can be considerably reduced after a few orienting runs. 



• ELIMINATES EIGENVALUE TRACKING 

• REQUIRES ONLY SINGLE ITERATION 

• REQUIRES SOLUTION OF ONLY ONE SET OF M ANO « 
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NUMERICAL METHODS FOR DIRECT SOLUTION 


The flutter Mach number and flutter frequency Mp and up may be solved for 
by Newton's method. The iterative scheme for Newton's method is 



where k is the iteration number and [Jr] is the Jacobian matrix given by 


u k ] = 


D 


kR,M t 


kl *M 


kR,co 


D 


kl 


FJ 


The Jacobian [J^l is expensive to compute because the evaluation of the aerody- 
namic matrix [AjJ is computationally intensive. Several quasi-Newton algo- 
rithms that approximate the Jacobian in various ways are available. 


NEWTON’S METHOD : 
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[J k ] EXPENSIVE TO COMPUTE 

QUASI-NEWTON METHODS: 

[Jk] APPROXIMATED IN VARIOUS WAYS 
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A BETTER QUASI -NEWTON METHOD 


A quasi-Newton algorithm is proposed that is more efficient than the algorithms 
currently available for determining the flutter Mach number and the flutter 
frequency. The numerical scheme is based on the hypothesis that approximating 
the derivatives of the flutter matrix [B k ] provides a more accurate approxima- 
tion to the Jacobian matrix [J k ] than directly approximating the derivatives of 
the characteristic determinant. The numerical scheme based on this hypothesis 
approaches Newton's method in its superior convergence characteristics with 
the same cost per iteration as the secant method. 

The derivatives [B k ] , m f and C B lJ *io F of the flutter matrix are approximated by 
following a reasoning similar to that employed in deriving Broyden's method 
(Johnston, 1982). Let AM k = M F ( k _i) - M F ( k ) and Aco k = w F ( k _i) - co F ( k ). The 
derivatives are approximated in the direction of the last move to satisfy 

t B k-ll = + t B k^*M F * A^k + t B kJ »co F * ^ u k 

and are assumed to be unchanged in the direction orthogonal to the last move. 


APPROACH 

• APPROXIMATE THE JACOBIAN INDIRECTLY BY APPROXIMATING 
THE DERIVATIVES OF THE FLUTTER MATRIX 

• UPDATE THE DERIVATIVES OF THE FLUTTER MATRIX ONLY IN THE 
DIRECTION OF THE LAST MOVE 

RESULT 

• A QUASI-NEWTON METHOD MORE LIKE NEWTON’S METHOD 
THAN OTHERS 

• THE FLUTTER MODE ALMOST A BYPRODUCT 


CD-88-31716 


1-330 


EFFICIENCY OF NUMERICAL PROCEDURE 


The direct solution procedure was demonstrated by performing flutter boundary 
calculations at various rotational speeds for two propfan rotor configurations. 
The first configuration consisted of eight identical blades. The second con- 
figuration was an alternately mistuned rotor with eight blades. 


The typical progress of iteration, for initial guesses for flutter Mach number 
and flutter frequency of 0.5 and 310 Hz, respectively, with the direct solution 
procedure and the conventional procedure, is shown in the first table. Recall 
that the conventional procedure relies on user interaction and judgment. The 
progress of iteration shown for the conventional procedure is typical. The 
direct solution procedure, in addition to being suitable for automation, is 
also more efficient as evidenced by the considerably smaller number of analysis 
steps. Thus both the CPU time and the analyst's time are considerably reduced 
by using the direct solution procedure. 

The results show that a fair initial guess would converge to the "exact" flut- 
ter point after about 5 to 10 flutter matrix evaluations. The second table 
shows the CPU times on the Cray-XMP required to obtain the flutter boundary for 
good initial guesses and poor initial guesses. The CPU times for one flutter 
eigenvalue analysis at a given set of Mach number and assumed frequency are 
also shown for comparison. With a good initial guess the flutter Mach number 
and the flutter frequency can be obtained for two or three times the cost of a 
single eigenanalysis . The direct solution procedure is much less expensive in 
terms of CPU time as well as analyst’s time than the conventional procedure, 
although precise comparisons have not been made. 


PROGRESS OF ITERATION 

(5280 rpm; BLADE SETTING ANGLE AT 0.75 RADIUS, 61.6°; a f =225°) 


COUNT 

CONVENTIONAL PROCEDURE 

DIRECT SOLUTION PROCEDURE | 


FLUTTER 

FLUTTER 

FLUTTER 

FLUTTER 


MACH 

FREQUENCY, 

MACH 

FREQUENCY, 


NUMBER, 


NUMBER, 

w, 


M 

Hz 

M 

Hz 

1 

0.500 

310.0 

0.500 

310.0 

2 



267.5 

.499 

310.0 

3 



268.9 

.500 

313.0 

4 



268.9 

.701 

289.7 

5 

.700 

268.9 

.590 

287.7 

6 



299.8 

.641 

293.6 

7 



298.9 

.642 

294.1 

8 



298.9 

.641 

294.1 

9 

.616 

286.3 

a .641 

294.1 

10 



290.1 



11 



290.5 



12 



290.4 



13 

.640 

292.9 



14 

.640 

293.9 



15 

.640 

294.0 



16 

.641 

293.9 



17 

.641 

294.1 



18 

a .641 

294.1 




Converged. 
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CPU TIME FOR AUTOMATED FLUTTER ANALYSIS 



EIGENSOLUTION 
AT A SINGLE 
SET OF MACH 
NUMBER M 
AND 

FREQUENCY « 

DIRECT SOLUTION TO FIND 
M F AND <<?f 


POOR 

INITIAL 

GUESS 3 

CPU TIME, SEC 

TUNED ROTORb-8 BLADES; 

6 MODES/BLADE; 5280 rpm; <r r =225° 

4.332 

10.356 
(M=0.70; 
<o=310 Hz) 

22.146 
(M = 0.45; 
w = 340 Hz) 

Ml STUN ED ROTORc-8 BLADES; 4 GROUPS; 
TWO MODES/BLADE; 5190 rpm; a r =90° 

10.020 

22.084 
(M=0.65; 
w=310 Hz) 

31.970 
(M=0.5; 
(o=340 Hz) 


aiNITIAL GUESSES ARE GIVEN IN PARENTHESES. 
b“EXACT” M F = 0.641 AND “EXACT” a> F =294 Hz. 
C “EXACT” Mp=0.718 AND “EXACT” <o F =285 Hz. 
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I 


ACCURACY OF JACOBIAN 


So that the accuracy of the approximate Jacobian could be judged, the nonlinear 
equations were solved by Newton's method, by the present numerical procedure, 
and by alternative quasi-Newton methods such as the multipoint secant method 
(implemented in IMSL routine ZSCNT), the modified Powell algorithm (implemented 
in IMSL routine ZSPOW), and Broyden's method (Johnston, 1982). For the cases 
tested, the present procedure outperformed all three alternative methods in 
terms of efficiency. Even though the characteristic determinant is never 

calculated in the present procedure, the variation of the absolute value of 
Dfc with each iteration is shown in the graph so that the procedure can be com- 
pared with Newton's method and the multipoint secant method. For these cases, 
the "exact" flutter Mach number was 0.641, the "exact" flutter frequency was 
294 Hz, and the initial values for Mp and cop were 0.65 and 330 Hz, 
respectively. The determinant value has been scaled so that 1.0 < Dq < 10.0, 
where Dq is the characteristic determinant at the beginning of iteration. 

The iteration history for the current numerical procedure closely resembles 
that for Newton's method, indicating the accuracy of the approximation pro- 
posed here for the Jacobian matrix. In contrast, the secant approximation for 
the Jacobian matrix requires almost double the number of iterations. 


5280 rpm; BLADE SETTING ANGLE AT 0.75 RADIUS, 61.6°; <r r =225° 

□ NEWTON’S METHOD 
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RANGE OF CONVERGENCE 


The range of convergence is an important factor in any iterative procedure 
since it has an important effect on how closely the initial solution must 
approximate the final solution. The graph on the left shows the number of 
flutter matrix evaluations required for convergence; the initial guesses for 
Mp are varied and the initial guess for cop is fixed at 310 Hz. The range of 
Mach number convergence is from 0.2 to 0.8. The graph on the right similarly 
shows the number of flutter matrix evaluations required for convergence; the 
initial guess for Mp is fixed at 0.65 and the initial guesses for up are 
varied. The frequency range of convergence with the direct solution procedure, 
230 to 350 Hz, is slightly larger than that with the secant method. From 
these graphs it can be stated that the present procedure has a large range of 
convergence. 


□ SECANT METHOD 
A DIRECT SOLUTION PROCEDURE 



ASSUMED FLUTTER ASSUMED FLUTTER 

MACH NUMBER, M F FREQUENCY, co F , Hz 
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SUMMARY 


A direct solution of the equations of motion is demonstrated to be a reliable 
automated flutter analysis procedure if steady aerodynamic deformations are 
ignored. The direct solution procedure replaces the inner-outer iteration loop 
of the conventional procedure by a single iteration loop. A numerical proce- 
dure, based on an accurate and efficient approximation to the Jacobian matrix, 
is presented. The procedure is straightforward in concept, and results for 
test cases show good convergence properties. Since the procedure is iterative, 
it is particularly suitable for design optimization. As the optimal design is 
evolved, the flutter solution is expected to change incrementally from design 
to design, so that the previous solution provides good estimates for the cur- 
rent solution. 


• DEVELOPED A QUASI-NEWTON METHOD FOR 
DETERMINANT ITERATION 

• AUTOMATED THE PROPFAN FLUTTER ANALYSIS BY 
DIRECT SOLUTION 

• DEMONSTRATED GOOD CONVERGENCE AND EFFICIENCY OF 
DIRECT SOLUTION METHOD WITH ADVANCED AERODYNAMIC 
MODEL 
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